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ABSTRACT 

In order to study multj-dimensional unstable detonation waves, we have developed a high 
order numerical scheme suitable for calculating the detailed transverse wave structures of 
multidimensional detonation waves. The numerical algorithm uses a multi-domain approach 
so different numerical techniques can be applied for different components of detonation waves. 
The detonation waves are assumed to undergo an irreversible, unimolecular reaction A — > 
B. Several cases of unstable two dimensional detonation waves are simulated and detailed 
transverse wave interactions are documented. The numerical results show the importance of 
resolving the detonation front without excessive numerical viscosity in order to obtain the 
correct cellular patterns. 
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Introduction 

Detonation waves are intrinsically multi-dimensional unstable phenomenon as demon- 
strated evidently by the experiments of Oppenheim [1] and White [2], Since then, the simple 
steady Chapman- Jouquet theory [3], [4] has been re-examined for its limited explaination 
of most of the multi-dimensional features seen in real world detonation waves. The unique 
characteristic of multi-dimensional detonation waves are their cellular patterns which are 
the trajectories recorded on the wall by the transverse waves structures. Those transverse 
wave structures consist of so-called “triple points” which are of three shock configurations 
(an incident shock, a reflected shock and a Mach stem plus a contact discontinuity) [5]. The 
formation of those transverse wave structures moving along the main precursor detonation 
front attracted much attention. These phenomena raised the interest of experimentalists 
trying to measure the cell size of the cellular pattern [6]. 

There are many aspects in the investigations of detonation waves. Most important is 
the formation of a compressible detonation wave from laminar deflagration waves [7], [8] 
- a process called the DDT (Deflagration to Detonation Transition) process. In a DDT 
process, turbulent boundary layers are recognized as playing an important role in flame 
acceleration and the buildup of compressible pressure waves in front of flame fronts. The 
latter eventually will cause ‘explosions in explosions’ in the reacting flows and generate 
detonation waves. Applied mathematicians have made different attempts to identify the 
mechanism in the formation of triple points in the context of simplified mathematical models. 
In [9]-[l 1] Erpenbeck first used normal mode analysis on the linearized Euler system to study 
the stability of multi-dimensional detonation waves. Later, in [12], Strehlow introduced the 
concept of acoustic ray trapping to study the formation of Mach stems; this idea has been 
generalized by Majda [13] using high frequency asymptotics. With the rapid advent of 
modern computing capability, another available avenue in studying detonation phenomenon 
is by direct numerical simulations. 

In this paper, we will present a high order hybrid method in order to study two dimen- 
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sional detonation waves, which will be able to resolve the detailed transverse wave structures 
of multi-dimensional detonation waves. The numerical simulations done in the paper are for 
an idealized model of chemical reaction in which the reactant species is in an irreversible, 
unimolecular reaction A — > B with finite Arrhenius reaction rate. It is evident that this 
model can not represent all the effects of realistic chemical kinetics on the cellular structures 
observed in lab experiments. However, an accurate numerical solution of this simplified 
model will provide a better understanding of the physics involved in the onset and evolution 
of detonation waves and a verification of current mathematical theories on detonation waves. 

There are three basic characteristics of detonation waves: (1) a strong precursor detona- 
tion front; (2) Mach stem configuration of the “triple points” and transverse wave structures; 
(3) stiff chemical reactions. The flow field can be divided into regions of highly irregular and 
steep gradients near the detonation front and regions of strong but smooth pressure waves. 
Near the shock fronts, strong vorticity fields are expected from the roll-up of slip lines. The 
temporal changes of thermodynamic and chemical compositions also vary dramatically from 
region to region. We will construct our numerical schemes according to these characteris- 
tics of multi-dimensional detonation waves, therefore , it is not surprising that the resulting 
numerical scheme is of a hybrid type. 

The reaction rate depends on the flow temperature exponentially through the Arrhenius 
relation. Accurate computations of the flow field are extremely important in producing the 
correct chemical reactions and thus the correct cellular structures. Traditional shock captur- 
ing schemes, designed to smooth shock and contact discontinuities, introduce a considerable 
amount of numerical viscosity near those discontinuities. They have been shown to have 
a tendency to distort the real chemical reaction processes. In [14], the widely used P.P.M. 
high-order Godunov scheme was found to produce nonphysical weak detonations. Also in 
[15], the ENO finite difference scheme was shown to yield wrong detonation speeds in one 
dimensional ZND simulations. All these facts point out the importance of designing numer- 
ical methods without excessive numerical viscosity. Among the attemts of simulating two 
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dimensional detonation waves numerically, Oran et al has done a series of simulations with 
FCT schemes and a phenomenological chemistry model was used to solve the stiffness in the 
system, see for instance [16]. Another early result was obtained by Taki and Fujiwara [17] 
where a quasi first order numerical scheme was used with a two front model suggested in 
[18] to model the chemistry. Most recently, an improved version of P.P.M., which uses the 
location of the detonation in evaluating the numerical fluxes, gave improved results in one- 
and two-dimensional detonation simulations [19]. 

The work reported in this paper is part of a larger research project which intends to 
understand different aspects of detonation waves, including the DDT process and the effects 
of turbulence and the detonability and detonation failures. In Section 1, we will introduce the 
governing equations for reacting flows and its formulation in general curvilinear coordinates. 
In Section 2, we discuss the hybrid numerical scheme using multi-domain approaches, and 
the ENO finite difference methods and Chebyshev collocation methods and shock tracking 
methods will be introduced. These numerical techniques will be used in the framework 
of multi-domain to fit the properties of 2-D detonation waves. In Section 3, we consider 
the treatment of interface conditions between different numerical schemes and boundary 
conditions and the smoothing techniques for the detonation front. In Section 4, we validate 
the hybrid numerical scheme and test the effects of smoothing of the detonation front on the 
cellular pattern of detonation waves. Then, we present the main results of this paper, the 
simulations of several cases of two dimensional detonation waves. Detailed analysis of the 
results will be discussed. Finally, in Section 5 we give a conclusion and the plan for future 
works. 
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1 Governing Equations 


Consider two dimensional detonation waves in an infinite channel moving from left to right 
into unreacted gas mixtures and the channel is denoted by f l, 


W W 

f) = (-00,00)*[-y, y] 


( 1 ) 


where W is the channel width. 

The governing equations for reacting detonation waves with s species and p reaction steps 
are the following Euler equations, 

du df{u) dg(u) 


dT + dx 


+ 


dy 


= *(u) 


( 2 ) 


where 


U = (p,pu,pv,pe,p u ---,p,) T 


f( 


u) = f/w, 


pu 2 + p, puv, pu(e + -) 


g(u) = ypv,pvu,pv 2 + p,pv(e + 
$(u) = (0, 0, 0, 0,u>i, • • • , w s ) T 


,Piu, -,p s uj 
-),Piv,---,p s v j 


where (u, v ), p, p, e are velocity vector, density, pressure, and total specific internal energy, p, 
is the mass density of the ith species. = pMiYTi= i ^ ere ^ ' s molecular weight 

of species -i and r'.j is the stoichiometric coefficient for the j — species in the i-th reaction 
step, and ^ denotes the change rate of the ith reaction progress variable which is assumed 
to obey Arrhenius’ rule. In the case of one-step A — ► B irreversible reaction, i.e. s = p = 1 
and 

w = —KpX exp(— ^-) (3) 


where A = is the mass fraction of reactant. If we assume exthermoic reaction, the specific 
internal energy 


e 


( 7-1 )p 2 


it 2 + v 2 

H “ I - ^ Q 


(4) 
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where Q is the specific heat formation and 7 = 1.2 is the ratio of specific heats. 

All quantities above have been non-dimensionalized by the initial states in the unreacted 
gas mixture in front of the detonation fronts. They are given as follows (“*— indicates 
nondimensionalization and “0” subscript denotes states of unreacted gas ) 


P 

P 

u 

T 

E+ 

Q 

t 


_E_ 
P 0 

P 0 
T_ 

A 

RT 0 

_2_ 

RT b 


Where t* = ;1V , l* = x\ is the half reaction distance which a mass particle will travel from 

V 2 

the detonation front before half-depletion of the reactant occurs. 

Consider a general curvilinear coordinates {£,r],t),£ = £(x,y,T),r] = y(x,y,T),t = 
t(x,y,T), the governing equation (2) becomes 


dU <9F(U) 0G(U) 

~dt + d£ + dy 


= *(U) 


(5) 


where U = ^/gu, 


F(U) = 


U'y/gp \ 


( u 2 \/gp \ 

U'y/gpu + 


U 2 y/gpu - ytf 

U l \/gp v ~ x nP 


U 2 y/gpV + Xtf 

U'yjgpe + (y v u - x, t v)p 

G(U) = 

U' 2 y/gpe + {~ViU + x^v)p 

U'y/gpx 


U 2 ^g P 1 

U l y/gp s y 


^ U 2 y/gps / 


( 6 ) 


U) = (0, 0, 0, 0, y/gu>l, • • • , \/gu, s ) T 


where ^fg — x^y n — x v y $ is the Jacobian of the mesh transformation, I/ 1 = 6 + u tx + 
v£ y , and U 2 = rj t + ur Jx + vy y . The spatial discretization is done in the (£,77) space for 
equation (5) and the range for £,77 are both [-1,1]. 

In most of the numerical simulations done in the paper, for initial conditions we use 
the Chapman- Jouguet ZND steady solution which can be obtained by solving the Rankine- 
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Hugoit equation between the state in front of detonation front and all states behind the 
detonation front [5]. 

2 Hybrid Numerical Algorithm with Domain Decom- 
position Technique 

In this section we describe the hybrid numerical scheme for detonation waves. The compu- 
tational region is composed of the detonation front moving to the right and the rear piston 
boundary and upper and lower solid walls. It will be subdivided using multi-domain tech- 
nique with the detonation front as the right most boundary (see Figure 1). Three different 
numerical techniques will be applied in different parts of the computational region. They 
are 

• Shock tracking algorithm for the detonation front; 

• High Order ENO finite difference scheme in the subdomain which contains the reflected 
shocks and contact discontinuities along the detonation front; 

• Chebyshev collocation methods for the strong vorticity and pressure fields from the 
interaction of transverse waves along the detonation front. 

We give a brief descriptions of each of the three numerical techniques. 

2.1 Chebyshev collocation methods 

In the computation domain (£,7/) € [ — 1 , 1]X[— 1, 1], to approximate any quantity /(£, /;), we 
consider its Chebyshev collocation interpolant lNf(£,v) which is defined as (assuming that 
polynomials are of the same order N in both the £ and rj directions): 

'«/((,?) = E (?) 

where 4>j(£) = — — , Co = = 2 and cj = 1, for 1 < j < N — 1. Here £, = cos ^ 

are the Chebyshev-Gauss-Labotto points, and Tjv(£) = cos(A^ cos _1 (£)) is the Chebyshev 
polynomial of the first kind. 
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The derivatives of /(£,*}) can be approximated by those of be. 

= E ( 8 ) 

0<i,j<N 

= E Mm)[ E 

0 <j<N 0 <i<S 

where the inner summation in the square bracket can be evaluated either by matrix-vector- 
multiplication on a vectorized machine or by Fast Fourier Transformation. The latter method 
only involves 0(N 2 log N) operations for the computation of all »?/)• Similar procedures 
can be obtained for f v (tk,Vi)- 


2.2 High Order ENO Finite Difference Methods 

To apply the ENO finite difference method to (5), the spatial derivative is discretized by 
conservative numerical flux differences: 


| , ^*J+2 _ iTi l IT. .) 

~ir + a ? + a , 


( 9 ) 


where Uij is the numerical approximation of (5) using the method of lines. In order to com- 
pute the numerical fluxes F i+ Lj,G i j+ 1 , first the primitive functions for F, G are approx- 
imated by piecewise polynomials using the ENO adaptive stencil idea of Harten [20], and 
then the derivatives of those polynomials are evaluated at the edge-centered mesh points 
(^ + l, T}j), (&, Vj+i) to produce the numerical fluxes. For technique details on the construc- 
tion of ENO fluxes to the system of equations (5), we refer the reader to [20], [21]. 

We need the eigenvalues and eigenvector on the Jacobian of fluxes A = gjj,B = aTj 
for the characteristic decompositions in order to define ENO numerical fluxes. If ^ = 
4- + u-§~ + v4~ denote the material derivative, then we have the following eigenvalues and 

dt ' ox dy 1 

eigenvectors, 
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Eigenvalues 
of A 

Eigenvalues 
of B 

Left Eigenvector 
of A 

Right Eigenvector 
of B 
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Hi 
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Dt 
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Hi 

Dt 

Du 

Dt 
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Dt 

Dr 1 
Dt 
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1 4 T 

T-'r 4 


and 


where a = 


1, = 

(l,-n x/ oa,-n y /9<z,0,0) 

ri = 

1 ( | n x Tly 1 

2 ' ’ pa 5 pa ’ a 2 ’ 

0) T 

12 — 

(0, — n y ,n x , 0,0) 

r 2 = 

(0, — rc y , n x , 0, 0) T 


I3 = 

(1,0,0, -<z 2 ,0) 

r 3 = 

p 

O 

P 

"l 

y- 

0 

H 

(10) 

1 4 — 

(0,0, 0,0,1) 

r 4 = 

(0,0,0,0,1) T 

i s = 

(l,n x pa,n y pa, 0 , 0 ) 

r 5 = 

1(1, 

2 ' ’ pa 7 pa ’ ’ / 


\/? is 

the sound speed and n x 

— — £* — n Zy 

\Al +«? ' y \Al+^y 

for A and n x = 


^4+u* ’ Uy = and the transform matrix T = jg with v = (p,u,v,p, A) being the 

primitive variable, 


1 2 i ( u2 + y2 ) -(7 _1 ) u —(7 — l)f (7-1) —(7 — 1 )Q \ 

-} i 0 0 0 

T= -f 0 i 0 0 

1 0000 

i 0 001/ 

2.3 Tracking Algorithm for Precursor Detonation Front 


( 11 ) 


The detonation front will be represented by a continuous curve 


x = 


-IT 

~2~ 


1T 

< y < — , * > 0 


( 12 ) 
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and the normal of the front denoted by n = ( n x ,n y ) points to the unreacted gas, which can 
be computed by 


n r = 


n y = 


\/l + ( x y)‘ 


(13) 


-x v 


v/ 1 + (^) 2 

By the Hygen’s principal, the detonation front will propagate in its normal direction with 
speed 


D n = 


x t 


(14) 


\A + ( x v) 2 

Following the procedure proposed in [23] [22], we can write the Rankine-Hognoit condition 
across the shock front as follows 


Po{u n ,0 - D n ) = pi(u ni i-D n ) 

Po( u n,o — D n ) 2 + po = pi(u n ,i ~ D, i) 2 + pi 
^T^ + 2 ( u n,o ~ D n ) 2 + XoQ = + !(«„,! - D n ) 2 + X\Q 


(15) 


PoAo(^n,0 I^n) — /^lAi(u nj i ^n) 
where “0” denotes the state in front of detonation front and “1” the states behind the front. 
And u n = (u,v) o n is the normal velocity on the front. Equation (15) relates the states 
in front of the shock and behind the shock. To close the system, we need to impose the 
continuity of tangential velocity, i.e. 

u t ,o = u M . (16) 

where u t = (— n y ,n x ) o (u,u) is the tangential velocity. 

From (15) (16), we can solve the quantities (p lt p\,u Hi i, u t) i, Ai) behind the shock in terms 
of those in front of the shock (po> Po ? u n,o> ^i,o, Ao). Using the notation in [9], we have the 
following 

(7+ l)/c 2 


P i — Po 


P 1 = Po 


( 7 k 2 + 1)(1 - in) 

( 7 K 2 + l)(l + W'f) 

7 + 1 


(17) 


«n,l = Un, 0 - («n.O ~ A i)(l - — ) 

P\ 
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u t,\ — 0 

Ai — A 0 


where w = \\ and k = I Pa tin i 2 i j s the Mach number of the shock front relative to the 

1*1 Co 

unreacted gas and c 0 = \f^~ “ y/l ls the sound speed in the unreacted mixture. 

In order to derive a time evolution equation for the shock front, we define the time 
differentiation along the shock front j t — Jj + x t (y,t)j for any fixed y € [-W/ 2 , W/ 2 ]. By 
applying ^ on both sides of ( 17 ) and separating the terms which involves x ti (y, t ), we obtain 
the following 


dpi 

dt 

dpi 

dt 

du i 

dt 

dv\ 

dt 

dXi 

dt 


— c\Xtt + d\ 

= C 2 x tt + ^2 

— c 3 >u x tt + dz }U 

— ^3 }V Xtt + d^y 
= C 4 Xtt + J 4 


( 18 ) 


where 


Cl = 

di = 

c 2 = 

d 2 = 

C 3 ,u “ 
= 


^ 3 ,v 

d 3,« 


(l+ 7 )(l+(^) 2 ) 

<*Po 2 7 k 2 + I-7 2^o ,, M , 2d 

(T+t)(i + (^) [2tA ' + ‘ 5,1 


dt 


1 + 7 
4(7 + l)/c 2 
r[(7 - l)« 2 + 2] 2 

P\ dpo 2(7 + l)/>o« 2 f 27 V f _ & 
/?o dt ' [(7 — 1)k 2 ■+■ 2] 2 r S l 

C3,n x y^3,t 

1 + K) 2 

^3,n + Xydz t i 

1 + (Xy) 2 

3,n + C 3 it 


1 + (x v ) 

_ ~Xy^3, n + t?3,t 

1 + (Xy) 2 
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93, t — T t — Xy t U\ 


<*." = ^rj( 1 + ^ c «-° 

_ (y^jy-2 2 5, 

d 3 ,n - X yt V!+ (7 + 1)k2 N t ( 7+1 ) r ’ 


U = 0 , g 4 = 


d A 0 

dt 


and S = (1 + (x y ) 2 )cl, N = u 0 - x y v 0 , T = x y u 0 + v 0 and N t = T t = S t - jjj. 

Using = |X^L. = T^, equation (18) can be rewritten for conservative variable u, 

dui 

m x 


dt 


— Tcx tt + T d 


(19) 


where c = (cj, c 2 , c 3iU , c 3|V , c 4 ) T , d = (di, d 2 , d 3 , u , d 3) „, d 4 ) T . 

Now on a fixed point on the shock front we consider the characteristic component of 
equation (2) in the normal direction n = ( ps ^ 


5 u d ~ d . 

at an at g 


( 20 ) 


where t is the tangential direction along the shock front, f = n x i + n y g, g — — n y f + n x g. 


The eigenvalues and eigenvector of Jacobian matrix A = are 


U n C, 1 i n , U n , u n c 


i a r,i 2 r,. ..,i 5 r 


(21) 


( 22 ) 


and i/d* Tiyi) . 


where !,-,* = 1, ■ • • ,5 are given in (10) with n x = ~^—,n y - ^> 1+( ^ )2 

Along the normal direction of the shock front, the (u n + ^-characteristic field approaches 
the shock front from behind, therefore a compatibility condition can be obtained by consid- 
ering the characteristic combination of equation (19) 

,duj 


hr- 


dt 


lsTcxt* + UTd, 


(23) 


thus yielding the following ODE for the shock speed 


x tt = H(u u \io,x y ,x t ,Xt y ) 


(24) 
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where 


ur^-urd 

I 5 Tc 

iu,_ u d 1 1 dV x 

dt ~ 1 dt v ^ + v /g~aT 

Qp/TJ\ 

— $ 1} + ^(U) is the residual computed in (£, tj) space behind the shock 

3 Interface Conditions and Smoothing of Detonation 
Front 

3.1 Interfacial Conditions between Subdomains and Boundary 
Conditions 

We first discuss the issue of interface and boundary conditions for the hybrid numerical 
scheme. A correct interface coupling between the ENO finite difference method and spec- 
tral method is crucial for the global stability and accuracy of the scheme. The basic idea 
behind the treatment of interface and boundary condition is the propagation of information 
along characteristics of the hyperbolic systems [24]. We will consider the following situation 
separately. 

Case 1. Interface between subdomains 

On a typical interface between two subdomains, say T between fit and f2 r in Figure la, 
there are two types of points, i.e. interior point I and cross point T. 

(a) Interior points I 

Let v = (P,u,v, p, A) T be the primitive flow variable, v^v 7 " denote the solutions com- 
puted for the time step C +1 in fy, fl r . Sr,i denotes the normal speed of the interface at point 
I with the normal direction n = (n x ,n y ). The eigenvalues and eigenvectors of the Jacobian 
matrix A are given in (10) and (22), and the corresponding characteristic variables are 

w i = p - apu n (25) 


and 


and ^ = 
front. 
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w 2 

= Ut 

(26) 

w 3 

= p-a 2 p 

(27) 

w 4 

= A 

(28) 

Ws 

= p + apu n 

(29) 


where denotes an average state between v ; and v r , for instance the Roe-average [28]. 

In order to update v at point I for the time step t n+1 , we make the following correction 
on Wi, 1 < i < 5 based on the sign of the difference between the eigenvalues and the normal 
speed of the interface Srj, i.e. 

corrected j ^ ^ i Sr, I < 0 

* -\ w\ if A, - S r , i > 0. 

Finally, we set v r = v ( = y C0TTected = corrected 

Remark. In the case that the interface is the shock front, v l should be computed by the 
Rankine-Hogoniot conditions (17). 

(b) Cross points between the wall boundary and interior interface - “T” points 

In order to update the solution for point “T” in Figure la we have to consider the 
information which comes from both subdomains ft/ and ft r and also the role of the wall. 
Characteristic surfaces approaching the point “T” from several direction can be used to 
form a closed system to determine this solution [25] [26]. Thus, such an approach is not 
unique and based largely on the experience. 

First let us consider the characteristic surface which is tangential to the interface T with 
normal n = (n r ,n y ) = (1,0). The corrected characteristic quantities w\,w 4 ,ws can be 
obtained as follows, assuming Ai,A 4 ,A 5 all positive (otherwise, replace subscript “1” by “r” 
for each negative values of A, — Sr,T, * = 1,4,5), 

w-y = w } = p — apu n (Jl) 

w 4 = W 4 = A (.} 2) 

tof = = v' + ml (33) 
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where subscript “x” indicates the normal direction of the characteristic surface and u n = uon. 

On the other hand, we know that the entropy s remains constant along the characteristic 
direction corresponding to the eigenvalue u n = u o n, we can correct the entropy s as follows 


_ corrected 


— | ^ Un ~ $ r ' T > ^ 

| otherwise 


Next, we consider the characteristic surface approaching the wall at point “T” with 
normal n = ( n x ,n y ) = (0, 1) - (top wall) or (0, —1) - (lower wall). On the top wall, u n — 
v = 0, so A/ = u n — a — ~a and we have the first characteristic field approaching the wall 
from the computational domain. Therefore, we can correct the first characteristic variable 
w\ using the results from either fi/ or Cl T , i.e. 


^ corrected ( P T 'f $r,T ^ 0 

1 1 [ p T T apu T n otherwise 


Finally, we solve for all four primitive quantities of the flow as follows 

v = 0 (36) 

A = w 4 (37) 

p = (w* + w% + 2iUj)/4 (38) 

P = (f)’. (39) 

The treatment of the cross point on the lower wall can be done similarly. 

Case 2. Wall Boundary Conditions 

For a point on the wall which only belongs to one subdomain such as point “B” in Figure 
la, n = ( n x ,n y ) = (0, 1) or (0, —1) for top and lower walls, respectively. So u n = u o n = 
v = 0, therefore the wall itself is a characteristic surface on which no normal flow condition 
is enforced, i.e. v — 0. Additionally, for the top wall (lower wall is treated similarly) we have 
the following compatibility equation according to its corresponding outflow characteristic 
field, 


w 2 = —n y u 
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W 3 

= p — a 2 p 

(41) 

W 4 

= A 

(42) 

W5 

= p + apu n . 

(43) 


So, the solution at point “B” can be obtained as follows, 

v = 0,u = ——,P = u> 5 , A = w 4 ,p = (p l - w 3 )/a 2 . (44) 


3.2 Smoothing Technique of Detonation Front 

In order to evolve the detonation front, we need to compute the time derivative of its normal 
speed D n , which depends on x tt in (24). Thus, the accuracy of the detonation front depends 
on the residual of numerical solution at the detonation front. Numerical experiments show 
that the front will develop high frequency numerical instability if no smoothing is applied on 
the detonation front. In this paper, we test three types of smoothing on the detonation front 
and compare the effects of different smoothing on the cellular pattern of detonation waves. 
Because detonation waves are unstable in most case in the high frequency range, this is a 
dilemma for our numerical simulation. On one hand, we try to admit as many frequencies in 
the numerical solution as possible in order to obtain enough nonlinear interaction between 
different unstable frequencies; on the other hand, to maintain numerical stability some cut- 
off in high frequencies is needed for long time integration. So, the best that we can hope for 
is to obtain as fine a resolution as possible in our numerical simulations while maintaining 
the numerical stability. 

Smoothing Technique One - Averaging Solution on the Detonation Front 

The simplest way to eliminate high frequency on the shock front is to smooth the solution 
Ui in //(uj, u 0 , x y , x t , x ty ) on the right hand side of (24). 

Smoothing I. In equation (24) replacing Uj tJ by the simple average of neighboring solutions, 
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which of course will reduce the accuracy of the solution on the detonation front, 

.. . UiJ-i + 2uj j + u lii+1 

u i j * ^ • 


( 45 ) 


It can be seen that this averaging procedure will reduce the accuracy of the shock front to 
only first order, and have a larger dissipation than the other smoothing techniques suggested 
below. 

Smoothing II - Derivative Smoothing 

One of the most used ipethods in shock capturing scheme to control oscillations derivative 
limiting. The idea was first introduced in [27] to construct monotonicity preserving cubic 
spline. 

Let (x,, yi),i = 1 , ■ • • , n be n- discrete data points and xj < X 2 • • ■ < x n . The cubic spline 
n 3 (x) is a piecewise cubic polynomial which has a continuous derivative at nodes x, and 
satisfies the following conditions, 


n 3(^) = Vi 

ny*,') =y\ 


and 


/ ^i-i 4" A, - _ i 

Vi = 


A,_i + Ax, 


1 < i < n. 


2<i<n-l 


(46) 


where s, = Ax, = x, +1 — x,. The derivative at the end points is given by 


, (2Aj + A2 )^i — AjS2 
Jfi = - 


Ax! + Ax 2 

/ (2A„_1 4" A„_2)s n -I A n _jS n _2 

Ax„_i 4- Ax n _2 

However, the cubic spline defined as above is usually oscillatory and the monotonicity 
of the original data set will be lost. It is suggested in [27] that the monotonicity can be 
preserved except at a local extrema point by limiting the derivatives y" s, 


, , _ j min(max(0,t/-),3min(s,_ 1 ,s,)) if y\ > 0 

| max ( mz ' n (Q, j/,-), — 3mm(s,_i , s,)) if y[ < 0. 


( 47 ) 


So we suggest the following smoothing procedure for the detonation front, 
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SMOOTHING II. In equation (24) change x' y in H as follows, 

(*' y )i ♦ — (*!)/,•«,.•, 1 <i<n. (48) 

Smoothing Technique III - High Frequency Spectral Space Smoothing 

The third smoothing technique tested in this paper applies high frequency cut-off func- 
tions in the Fourier transform space of the detonation front. Assume that the discrete data 
Xi, 1 < i < n are equally distributed and data (x,-, y,) has been decomposed into Fourier 
modes, 

Vi = E yk e ' kx ' l <i <n (49) 

where y k = J £" =1 Vie~ ikxi . 

We multiply the Fourier coefficient y ^ by a decreasing factor < 7 * so that the high frequen- 
cies will be decreased. The modified j/ t is given by 

--i 

Vi «- (yi) filtered =: E <rky k e ikxi l < i < n. (50) 

*=-7 

Here we chose <r k so that it decays exponentially in terms of the wave number, 

<7fc = for \k\ < (51) 

where the constant /r is chosen so that is the machine zero and 2£ is called the order of 
the exponential filtering. 

4 Cellular Structures of Two Dimensional Detonation 
Waves 

4.1 Linear Stability Analysis and 2 D Detonation Waves 

In [9]-[l 1] Erpenbeck first used linearized normal mode analysis to study the stability of two 
dimensional detonation wave. With complex analysis technique, he analyzed the unstable 
modes of linearized Euler equations with respect to the steady state solution of plane ZND 


17 



detonation waves. The stability of the detonation front xj>(y,z,t), thus the whole flow field, 
is determined by the existence of poles of its time Laplace transformation 7?(r,e) for any 
transverse frequency f = a 2 + /? 2 . It is found that, for large wave numbers (high frequencies), 
the stability of ZND detonation waves depends on the quantity c[](x) — u 2 (x) where x is the 
distance measured away from the ZND steady shock front and co(x) is the frozen sound 
speed at location x and u(x) is the flow velocity there. Detonation waves with an one step 
irreversible A — > B reaction have been categorized [11] in terms of their short wave stabilities: 

Type D Cg(x) — u 2 (x) is decreasing as x moves away from the detonation front; 

Type I Cq(x) — u 2 (x) is increasing as x moves away from the detonation front; 

Type M Cq(x) — u 2 (x) achieves a maximum at some points x* between x = 0 and x = — oo. 

Type D has been shown to be stable for high frequencies, while Type I and M will be 
unstable for bands of large wavenumbers, e,y < e < e,>,i = 1,2, ••• where e = fpr and W 
is the channel width. The latter case means difficulty in attempting to simulate detonation 
waves through numerical calculations as it will never be possible to resolve all the unstable 
modes in the system with finite a number of mesh points. 

4.2 Validation of Numerical Scheme 

Computational Mesh Set-up and Initial Conditions 

The following notation will be used in all the computations: ndm - the number of sub- 
domains in the subdivision of the computation domains; - channel width; Isk - number 
of marker points on the shock front; (n,m) - size of mesh in a subdomain. In all the com- 
putations presented here, we take ndm = 9 and the total length of the channel to be 150^ 
with the detonation front as the right boundary of the solution domain. As the detonation 
propagates and curves, the interfaces between subdomains will also travel at a speed which 
is taken to be the averaged speed of the curved detonation front. 

In the first subdomain, we use third order ENO-LF schemes[2l] in order to resolve the 
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reflected shocks and contact discontinuities. A stretching function will be used in the x- 
direction so that the mesh will be clustered toward the detonation front. The stretching 
function is given by 


2 . sin 
C = 0(0 = -1 + - sin — 


-i 


*K + i) 


-i<0f <i 


where is the stretched mesh and the parameter c* determines the amount of stretching, 
we take a = 0.995. 

The right hand side of the first subdomain, being the detonation front, will be tracked by 
the track algorithm in Section 2.3. Appropriate smoothing will be applied on the detonation 
front for about every 20 iterations; the exact frequency depends on the strength of the deto- 
nation waves. Chebyshev collocation methods will be used in the remaining subdomains. In 
order to ease the CFL restriction from the nonuniform distribution of Chebyshev collocation 
points, another stretching function <fi(x,a) is used to produce a more uniformly distributed 
mesh in the Chebyshev collocation subdomains. We take 


r = *(o = 


sin 1 a£ 
sin -1 a 


where again £* is the stretched mesh and a = 0.999 so that Chebyshev-Gauss-Labotto points 
in £ space will be mapped to (* space with more uniform distribution. 

Finer meshes will be used for those subdomains closer to the denotation front. A typical 
mesh set-up is given in Figure lb (only the first seven subdomains are shown). 

The computation starts with the ZND steady state solution of plane detonation wave. 
To induce the transverse waves, we introduce perturbations either in the detonation front or 
in the flow variables themselves, or both. The type of perturbations used will be sine-like 
wave 

x(j/,0) =: x(y,0) + esin(7ry 2 ), —W/2<y<W/2 (52) 


or random perturbations (eran/()). 

Reflective solid boundary condition will be used in all computations. 


19 


Effect of Smoothing of the Detonation Front on the Cellular Patterns 

The procedure of smoothing on the detonation front basically introduces extra numerical 
viscosity in the whole scheme, therefore the smaller this extra viscosity is, the more reliable 
the numerical simulations should be. This argument can be extended to any numerical 
simulation of detonation waves. To evaluate the effects of smoothing of the detonation front 
on the cellular structure, we consider the ZND detonation wave with heat formation Q = 50, 
activation energy E — 50, and the overdrive / = 1.6. The channel width is taken to be 
W = 20/*. As we are not. interested in the detailed structure of the flow field, we take a less 
fine mesh - £° =1 (n,m) = (50, 50) + (24, 50) + (20, 50) + (20, 50) + (20, 40) + (10, 20) + (10, 10) + 
(10, 10) + (10, 10). The number of marker points on the detonation front is Isk = 160. The 
initial size of the subdomains will be 5/* x W, 5/* x W, 5/* x W, 5/* x W, 10/* x W, 15/* x W, 
25/* x W, 40/* x W, 40/* x W. Three tests are done to see the effects of smoothing on the 
cellular pattern. In all the numerical results reported in this paper, the detonation front has 
traveled at least 20 channel width for the cellular patterns reported. 

Test One Smooth I, II, and III 

In the first test we activate all three types of smoothing on the detonation front. Thus, 
strong numerical viscosity will be produced to stablize the front. But, keep in mind, even in 
this situation, we still track the detonation front and no difference across the front is used 
in the scheme. In Figure 2a, we record the pressure on the detonation front for time t — 
10f* — 20f*. A very regular and symmetric two cell structure is produced by the interaction 
of four different triple points. This correspond to a cell width 10/* - half the channel width. 

Test Two. Smooth II and III 

In the second test, we deactivated Smooth I which produces the largest numerical viscosity 
among all three types of smoothing. In this case, only one cell is present in the cellular pattern 
(Fig 2b) which corresponds to a cell width 20/*. 

Test Three Smooth III only 
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In the third test, we use only Smoothing III which uses high frequency attenuation in 
the Fourier transformation space for the shock front. The order of exponential cut-off in 
(51) is 12, thus yielding very slight numerical viscosity on the detonation front. The cellular 
pattern of the detonation front (Fig 2c) consists of two quite irregular cells, a larger cell 
and a smaller cell. Two major triple points dominate the cellular structure along with two 
secondary triple points. A close-up picture of the cell pattern is given in Fig 6a where we 
will further examine this case in more detail. 

These tests show us the sensitivity of the cellular patterns of the detonation waves to the 
amount of numerical viscosity in the scheme. Consequently, we have to be very careful in 
applying the right kind of algorithm for the simulations if we want to compare the numerical 
results with either theoretical predictions or experiment results. Presumably, less numerical 
viscosity will produce more reliable cellular patterns. From this point on, all numerical results 
presented will use only the Smooth III procedure (about every 20 iterations), which has the 
smallest amount numerical viscosity, in order to stabilize the evolution of the detonation 
front. There have been several situations which demonstrate that such smoothing is necessary 
or the computation will abort prematurely. 

Accuracy of Time Integration and Mesh Convergence Studies 
1) Comparison of Time Discretizations 

The stiffness in a chemically reacting system poses a lot of difficulty for numerical sim- 
ulations. There are basically two issues to be considered when choosing a time integrator , 
one is accuracy and another is CPU time efficiency. For the one reaction system tested here, 
we could us either a time splitting method as in [19] [16] or just an explicit Runge-Kutta 
type method. For the splitting method, the evolution of the solution can be split into two 
steps, the first step being an Euler step where the governing equation is solved without the 
chemical reaction production terms; the second step involves only the reaction term with 
the temperature field frozen at the value of the previous Euler step. The second step can be 
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explicitly solved here because only one reaction is considered. However, if more complicated 
chemistry is involved, such splitting will not avoid the stiffness problem in the simulation. 

We compare the results of the splitting method (which is at most second order) and the 
third order Runge-kutta method with the same resolution and spatial discretizations in the 
spatial directions. The detonation parameters are Q = 50, E = 50, / = 3; the mesh sizes 
are X)i=i ( n > m ) — (100, 100) + (34, 70) + (20, 40) + (20, 40) + (20, 30) + (10, 20) + (10, 10) + 
(10, 10) + (10, 10); and the channel width W = 10**. The initial size of the subdomains will 
be 5** x W, 5 ** x W, 5 ** x W, 5** x W, 10** x W, 15** x W, 25** x W, 40 ** x W , 40** x W. In 
Figure 3, we plot the pressure along the center of the channel at time T = 34 1*. The solid 
line is the result obtained by the third order Runge-Kutta method while the dots (o) are 
the results by the splitting method. We can see both results agree fairly well in most parts 
of the domain except that the splitting method fails to resolve the dip in the pressure. In 
the numerical tests given later, the third order Runge-Kutta method will be used in all the 
computations. 

2) Mesh Convergence Studies 

We use the same detonation parameter as above but with three different meshes in the 
spatial direction, the third order Runge-Kutta method is used in both cases. Mesh A 
is 5Zi=i( n i m ) = (100,100) + (34,70) + (20,40) + (20,40) + (20,30) + (10,20) + (10,10) + 
(10, 10) + (10, 10). Mesh B is £? =1 (n,m) = (120, 150) + (34, 70) + (20, 40) + (20, 40) + (20, 30) + 
(10,20) + (10,10) + (10,10) + (10,10). Mesh B is £? =1 (n, rn) = (1.60, 20°) H- (34, 70) + (20, 40) + 
(20,40) -f (20, 30) + (10,20) + (10, 10) + (10, 10) + (10, 10). The initial size of the subdomains 
will be 8 1 x W, x W, 5t x W, 5t x W, 10T x W, 12 1 x W, 25** x W, 40** x W, 40** x W . 
So in the first subdomain, mesh A is about 10 points per ** and mesh B is about 15 points 
per ** and Mesh C is about 20 points per **. In Figure 4, we plot the pressure along the 
center of the channel ((-) Mesh C, (o) Mesh B, (+) Mesh A). Close agreement can be seen 
among the results for all meshes. In the rest of the computation, we will use at least the 
resolution of Mesh B, which is about 15 points per half reaction distance. 
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4.3 Numerical Simulation of 2-D detonation Waves 

We present four cases of detonation waves using our high order hybrid scheme, each belonging 
to one of the three types of detonations in terms of the short wave instability. 

CASE I Q = 10, E = 50,/ = 1.2, channel width W — 10r - One cell pattern 

This is a case of Type M detonation which is unstable for high frequencies. The mesh 
used in this case is E®=i (","*) = (HO, 140) + (20,40) + (20,40) + (20,30) + (20,30) + 
(10,20) + (10, 10) + (10, 10) + (10, 10). The number of “Marker Points” on the shock front 
is Isk = 300. The initial size of the subdomains will be 5 1 x W,5t x W, 5 1 x W,5t x W, 
lor x W, 15** X W, 25 1 x VP, 40r x W,40t x W. This mesh gives a resolution of 14 
points per half reaction distance in the ENO domain. As a result of the high accuracy of 
the Chebyshev collocation method, we find that very good resolution of the flow field can be 
obtained with far less points. We start the computation with the ZND steady state solution 
with a sine-like perturbation (52) on the shock front with e = 0.15. 

In Figure 5a, we record the pressure along the shock front for time t = 10<* — 20£*. Only 
one cell is produced in this run which corresponds to two triple points along the detonation 
front. In Figure 5b, we contour five snapshots of the pressure, temperature, vorticity and 
mass fraction at time T = 30.5f*,31.5t*,32.5t*,33.5t*,34.5f*. In Figure 5c, we sketch the 
interaction of the triple points which will be typical for all the other later cases. We can 
see the evolution of the reflected shock from the pressure field; the contact discontinuity 
can be best seen from the vorticity field. Because the contour lines hardly distinguish the 
exact position of contact lines, we can use the temperature field to locate the position of 
the curving contact lines. In the first time snapshot ( T = 30. 5f*) of both Figure 5b and 5c, 
we see two reflected shocks (RS) A,B waves moving toward the center of the channel and 
to collide. Following RS - A, B are two contact discontinuities C^,C+ respectively, where 
the signs indicate opposite circulation of these two contact lines which produce opposite 
vorticities. In snapshot two (T = 31. 5f*), reflected shock A, B have emerged from the 
interaction, exchanged directions, and are moving away from each other. Notice that in 
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the field of mass fraction, the layer of unreacted gas is much thicker behind the incident 
shock than the one behind the Mach Stem. The chemical reaction in the layer behind the 
incident shock waves provide the energy for further development of reflected shocks A, B. 
The contact discontinuities Cf, C+ are now detached from their triple point configurations, 
moving downstream, and their tips are rolled up. At the same time a new pair of contact 
lines Ci,Ci emerge behind the reflected shock waves B,A respectively. In snapshot three 
( T = 32. 5<*), RS - B, A are moving toward upper and lower wall respectively and are ready 
to be reflected away from the wall with the contact lines Cf, Cf following them. In snapshot 
three ( T = 33. 5t*), RS /?, A have been reflected away from the walls and the contact lines 

, Ci detached from the triple point configuration while an another new pair of contact 
lines CaiCf are created behind B and A respectively. In the last snapshot (T = 34. 5<*), 
RS - B, A are ready to collide again which finishes one cycle of the interaction of these two 
triple points. 

Case II Q = 50, E — 50, / = 1.6, channel width W = 20^*. Two cells pattern 

This is a type M detonation wave which is unstable for high frequencies, and we use the 
mesh £? =1 (n,m) = (50,250) + (34,70) + (20,40) + (20,30) + (20,30) + (10, 20) + (10, 10) + 
(10,10) + (10,10). The number of ‘Marker Points’ on the shock front is Isk = 300. Fig- 
ure 6a shows a two cell pattern produced by the trajectories of four triple points. There 
is a larger cell with width approximately 10£* (half the channel width) and a smaller 
one with width 5 1*. Figure 6b contains six snapshots of the detonation at time T = 
20.25f*, 21. 5f*, 22t*, 22. 5t*, 23<*, 23. 5t*. The times were chosen so that the interaction of 
the triple points can be shown clearly in the contour plots. A random perturbation with 
magnitude e = 0.3 is used to perturb the shock front at T = 0. Four triple points are pro- 
duced (two major ones and two secondary ones) and the cell pattern is given in Figure 6a. 
Referring to the six arrows and depicted shock fronts in the sketch Figure 6c which corre- 
sponds roughly the six time snapshots in Figure 6b . In the first time snapshot T = 20.25<*, 
in the lower middle part of the channel, two triple points C , D are moving away from each 


24 


other just after collision and A has been reflected from the upper wall and is about to collide 
with triple point B. In snapshot two (T = 21.5 f ), in the upper part of the channel, two 
smaller triple points A, B collide and separate and, in the middle part of the channel triple 
point C is moving upward to collide with triple point A. Near the lower wall, triple point 
D is reflecting away from the lower wall. In snapshot three (T = 22 t*), triple point B is 
approaching the upper wall, triple points A and C collide, and triple point D keeps moving 
up away from the lower wall. In snapshot four (T = 22.5 1*), triple point B reflects away 
from the upper wall and is about to collide with triple point C , while triple points A and 
D are about to collide with each other. In snapshot five ( T — 23 t*), in the upper part of 
the channel, triple points B and C collide and triple points A and D approach each other. 
Finally, in snapshot six ( T = 23.5 t*), triple points B and C finish the collision and exchange 
directions while triple points A and D collide. 

Also notice that in Snapshots 4 and 5, the two pressure waves from the reflected shocks 
intersect with each other before the collision of the triple points happens along the detonation 
front (Snapshot 6). Such interaction will cause sudden reaction of any unreacted gas in the 
interior region and produce so-called “explosions in explosions’. 

Case III Q = 50, E = 50,/ = 3, channel width W = 10 1. One cell pattern and chaotic 
flow fields. 

This case represents strong heat release, large overdrive, and belongs to type I which is 
again unstable for a range of high frequencies. We use a mesh X)f =1 (n,m) = (120, 150) + 
(34, 70) + (20, 40) + (20, 30) + (20, 30) + (10, 20) + (10, 10) + (10, 10) + (10, 10); and the number 
of ‘Marker points’ on the shock front is Isk = 300. The initial size of the subdomains will be 
5 t x W, 5 1 x W, 8 1 x W, 5 1 x W, 10 1 x W, 12 1 x W, T5t x W, 40t x W, 40f x W. Only two 
triple points are produced in this case; an one cell pattern of the detonation front is given 
in Figure 7(a). In Figure 7(b), we have six time snapshots of the pressure, temperature, 
vorticity and mass fraction. Large vorticities are generated behind the detonation front and 
the flow field becomes quite chaotic. There are only two triple points along the detonation 
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front which produce a one cell pattern for the cellular structure with cell size W = IOC (see 
Figure 7a). 

5 Conclusion 

We have developed a high order numerical scheme which is suitable for computing detailed 
transverse wave structures of two dimensional detonation waves. The numerical algorithm 
uses a multi-domain approach so that different numerical techniques can be applied for 
different components of detonation waves. The propagation of waves across the interfaces of 
subdomains have been very smooth and the order of accuracy of the whole numerical scheme 
is only limited by the accuracy of the time integrator. Tracking of the detonation front avoids 
differences across the detonation front, thus avoiding excessive numerical viscosity in shock 
capturing schemes. The high resolution of the Chebyshev collocation method enables us to 
use far less grid points in most of the solution domain and yields great savings in the total 
CPU cost. The potential for using higher ENO finite difference schemes in the subdomain 
which contains the reflected shocks and contact discontinuities can be further exploited. 

We have shown that the cellular pattern of the detonation waves is affected by the 
accuracy of the detonation front and the amount of numerical viscosity, especially the amount 
of viscosity involved in the time evolution of the detonation front by the numerical scheme. 
We believe that this point should be well taken in the further investigation of detonation 
waves in order to have meaningful comparisons with experiment results. 

We have studied several cases of detonation waves with specific ratio 7 = 1.2, from small 
heat release (Case I) to large release (Case II and III) and small overdrives (Case I) to 
large overdrives (Case II, III). The numerical results successfully reproduced the onset and 
evolution of the transverse wave structures. The contact lines within triple points create 
large vorticity fields behind the detonation front which will distort and interact with the 
detonation front. The contact discontinuities from the triple points after their collisions 
convect downstream and generate vorticity downstream. Further work will be done by using 
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more realistic chemistry models so that comparison with experiment results will be possible. 
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Figure 1 (a) Multidomain set-up for the hybrid numerical scheme for detonation waves; (b) 
Mesh structures for 2-D detonation waves 
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Figure 3 Pressures along the center of channel obtained by time splitting (0) and third 
order Runge-kutta (-). 



Figure 4 Mesh convergence studies: Pressure along the center of channel with three meshes 
(-) 20 points per (o) 15 points per (*, (+) 10 points per 
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Figure 5 Case I, Q = \0,E — 50,/ = 1.2, VC = 10 (a) one cell pattern, (b) five snapshots 
of pressure, temperature, vorticity, mass fraction of reactant (from top to bottom) at T = 
30.5r, 31.5f*, 32.5f*, 33.5C, 34.5C, (c) sketch of the interaction of the triple points. 
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Figure 6 Case II, Q = 50, E = 50,/ = 1.6, W = 20 (a) two cell pattern, (b) six snapshots 
of pressure, temperature, vorticity, mass fraction of reactant (from top to bottom) at T = 
20.25C, 21.5C, 22C, 22. 5C, 23C, 23. 5C. (c) Tracks of the triple points. 
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Figure 6c 
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Figure 7 Case III, Q = 50, E = 50, / = 3, W = 10 (a) one cell pattern, (b) six snapshots 
of pressure, temperature, vorticity, mass fraction of reactant (from top to bottom) at T = 
33C , 33.5 1\ 34 f, 34.5C , 35C, 35.5C. 
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